The purpose of this script is to evaluate the spatial and temporal overlap of the different surveys used in the SEM.

Inventory datasets:

Annual SEM: - Fish: FMWT, Bay Study, DJFMP Beach seine (ISS) - Zoop/nutrients: EMP - Benthic: EMP Benthic - Flows: Dayflow

Monthly SEM: - Fish: - Zoop/ nutrients: - Benthic: - Flows: Dayflow

Load libraries

# tidy data
library(tidyverse)

library(DT)
library(leaflet)
library(sf)
library(janitor)

Fish surveys

Query the CDFW Fall Midwater Trawl and Bay Study data sets from the LTMRdata package.

library(LTMRdata)

Surveys <- LTMRdata::fish(sources=c("FMWT", "Baystudy"),
             species=NULL,
             convert_lengths=TRUE,
             size_cutoff=40,
             zero_fill = TRUE,
             remove_unknown_lengths = TRUE,
             univariate = FALSE,
             quiet = FALSE)
Surveys
## # A tibble: 9,223,361 × 25
##    SampleID  Taxa  Length Count Length_NA_flag Source Station Latitude Longitude
##    <chr>     <chr>  <dbl> <dbl> <chr>          <chr>  <chr>      <dbl>     <dbl>
##  1 Bay Stud… Acan…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  2 Bay Stud… Acip…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  3 Bay Stud… Acip…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  4 Bay Stud… Allo…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  5 Bay Stud… Alop…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  6 Bay Stud… Alos…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  7 Bay Stud… Alos…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  8 Bay Stud… Amei…     NA     0 <NA>           Bay S… 101         37.5     -122.
##  9 Bay Stud… Amei…     NA     0 <NA>           Bay S… 101         37.5     -122.
## 10 Bay Stud… Amei…     NA     0 <NA>           Bay S… 101         37.5     -122.
## # … with 9,223,351 more rows, and 16 more variables: Date <dttm>,
## #   Datetime <dttm>, Survey <int>, Depth <dbl>, Method <chr>, Tide <chr>,
## #   Sal_surf <dbl>, Temp_surf <dbl>, Secchi <dbl>, Secchi_estimated <lgl>,
## #   Tow_volume <dbl>, Tow_direction <chr>, Cable_length <dbl>,
## #   Tow_duration <dbl>, Tow_area <dbl>, Notes_tow <chr>
Surveys$Month <- lubridate::month(Surveys$Date)
Surveys$Year <- lubridate::year(Surveys$Date)
str(Surveys)
## tibble [9,223,361 × 27] (S3: tbl_df/tbl/data.frame)
##  $ SampleID        : chr [1:9223361] "Bay Study 1" "Bay Study 1" "Bay Study 1" "Bay Study 1" ...
##  $ Taxa            : chr [1:9223361] "Acanthogobius flavimanus" "Acipenser medirostris" "Acipenser transmontanus" "Allosmerus elongatus" ...
##  $ Length          : num [1:9223361] NA NA NA NA NA NA NA NA NA NA ...
##  $ Count           : num [1:9223361] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Length_NA_flag  : chr [1:9223361] NA NA NA NA ...
##  $ Source          : chr [1:9223361] "Bay Study" "Bay Study" "Bay Study" "Bay Study" ...
##  $ Station         : chr [1:9223361] "101" "101" "101" "101" ...
##  $ Latitude        : num [1:9223361] 37.5 37.5 37.5 37.5 37.5 ...
##  $ Longitude       : num [1:9223361] -122 -122 -122 -122 -122 ...
##  $ Date            : POSIXct[1:9223361], format: "1980-01-23" "1980-01-23" ...
##  $ Datetime        : POSIXct[1:9223361], format: "1980-01-23 07:49:00" "1980-01-23 07:49:00" ...
##  $ Survey          : int [1:9223361] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Depth           : num [1:9223361] 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 12.2 ...
##  $ Method          : chr [1:9223361] "Midwater trawl" "Midwater trawl" "Midwater trawl" "Midwater trawl" ...
##  $ Tide            : chr [1:9223361] NA NA NA NA ...
##  $ Sal_surf        : num [1:9223361] 20.5 20.5 20.5 20.5 20.5 ...
##  $ Temp_surf       : num [1:9223361] 11.5 11.5 11.5 11.5 11.5 11.5 11.5 11.5 11.5 11.5 ...
##  $ Secchi          : num [1:9223361] 86 86 86 86 86 86 86 86 86 86 ...
##  $ Secchi_estimated: logi [1:9223361] NA NA NA NA NA NA ...
##  $ Tow_volume      : num [1:9223361] 5603 5603 5603 5603 5603 ...
##  $ Tow_direction   : chr [1:9223361] NA NA NA NA ...
##  $ Cable_length    : num [1:9223361] NA NA NA NA NA NA NA NA NA NA ...
##  $ Tow_duration    : num [1:9223361] NA NA NA NA NA NA NA NA NA NA ...
##  $ Tow_area        : num [1:9223361] NA NA NA NA NA NA NA NA NA NA ...
##  $ Notes_tow       : chr [1:9223361] NA NA NA NA ...
##  $ Month           : num [1:9223361] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Year            : num [1:9223361] 1980 1980 1980 1980 1980 1980 1980 1980 1980 1980 ...
detach(package:LTMRdata, unload=TRUE) # detach because it conflicts with discretewq

Assess sampling density across space and time

SamplingDensity <- Surveys %>%
  select(Source, Method, Station, Latitude, Longitude, Year, Month) %>% 
  distinct(Source, Method, Station, Latitude, Longitude, Year, Month) %>% 
  mutate(Sampling=1)

summary(SamplingDensity)
##     Source             Method            Station             Latitude    
##  Length:55488       Length:55488       Length:55488       Min.   :37.50  
##  Class :character   Class :character   Class :character   1st Qu.:37.89  
##  Mode  :character   Mode  :character   Mode  :character   Median :38.04  
##                                                           Mean   :37.98  
##                                                           3rd Qu.:38.07  
##                                                           Max.   :38.56  
##                                                           NA's   :4      
##    Longitude           Year          Month           Sampling
##  Min.   :-122.5   Min.   :1968   Min.   : 1.000   Min.   :1  
##  1st Qu.:-122.4   1st Qu.:1992   1st Qu.: 4.000   1st Qu.:1  
##  Median :-122.2   Median :2001   Median : 8.000   Median :1  
##  Mean   :-122.1   Mean   :2000   Mean   : 7.168   Mean   :1  
##  3rd Qu.:-121.8   3rd Qu.:2010   3rd Qu.:10.000   3rd Qu.:1  
##  Max.   :-121.4   Max.   :2020   Max.   :12.000   Max.   :1  
##  NA's   :4
# Bay Study
Baystudy <- SamplingDensity %>% 
  filter(Source=="Bay Study") %>% 
  group_by(Station, Longitude, Latitude, Year) %>% 
  summarize(Sampling=sum(Sampling)/24) # 2 gears x 12 months = 24 sampling points per year
Baystudy
## # A tibble: 1,962 × 5
## # Groups:   Station, Longitude, Latitude [56]
##    Station Longitude Latitude  Year Sampling
##    <chr>       <dbl>    <dbl> <dbl>    <dbl>
##  1 101         -122.     37.5  1980    0.917
##  2 101         -122.     37.5  1981    1    
##  3 101         -122.     37.5  1982    1    
##  4 101         -122.     37.5  1983    0.917
##  5 101         -122.     37.5  1984    0.958
##  6 101         -122.     37.5  1985    0.917
##  7 101         -122.     37.5  1986    0.917
##  8 101         -122.     37.5  1987    0.917
##  9 101         -122.     37.5  1988    0.917
## 10 101         -122.     37.5  1989    0.667
## # … with 1,952 more rows
ggplot(Baystudy, aes(Longitude, Latitude))+
  geom_point(aes(color=Station, size=Sampling), pch=21)+
  theme_bw()+
  facet_wrap(Year~.)+
  scale_color_viridis_d()

# FMWT
Fallmwt <- SamplingDensity %>% 
  filter(Source=="FMWT") %>% 
  group_by(Station, Longitude, Latitude, Year) %>% 
  summarize(Sampling=sum(Sampling)) 
Fallmwt
## # A tibble: 4,951 × 5
## # Groups:   Station, Longitude, Latitude [180]
##    Station Longitude Latitude  Year Sampling
##    <chr>       <dbl>    <dbl> <dbl>    <dbl>
##  1 070         -122.     38.2  1992        7
##  2 070         -122.     38.2  1993        6
##  3 070         -122.     38.2  1994        5
##  4 070         -122.     38.2  1995        1
##  5 071         -122.     38.2  1992        4
##  6 072         -122.     38.2  1991        3
##  7 072         -122.     38.2  1992        7
##  8 072         -122.     38.2  1993        6
##  9 072         -122.     38.2  1994        9
## 10 072         -122.     38.2  1995        7
## # … with 4,941 more rows
ggplot(Fallmwt, aes(Longitude, Latitude))+
  geom_point(aes(color=Station, size=Sampling), pch=21)+
  theme_bw()+
  facet_wrap(Year~.)+
  scale_color_viridis_d()

Zoop surveys

Query the EMP data from the Zooper package. Borrow code from “zoop.Rmd” by S.B.

library(zooper)

Load core stations

download.file("https://portal.edirepository.org/nis/dataviewer?packageid=edi.522.7&entityid=71dd301f30a2bc2e40f5da573dde9f97", destfile = file.path(tempdir(), "zoop_station_lookup.csv"))

zoop_stations<-read_csv(file.path(tempdir(), "zoop_station_lookup.csv"))%>%
  filter(Core%in%c(1, 2))

Load zoop data

zoop_data<-Zoopsynther(Data_type="Community", Sources="EMP", Time_consistency = TRUE)%>%
  filter(Station%in%zoop_stations$StationNZ)%>%
  mutate(Month=lubridate::month(Date))
## [1] "No disclaimers here! Enjoy the clean data!"
detach(package:zooper, unload=TRUE)

Remove stations that aren’t continuous

zoop_data_continuous<-filter(zoop_data, !Station%in%c("NZEZ6", "NZEZ2", "NZD16", "NZD06", "NZ080", "NZ042"))

# Make sure the right number of stations was removed
length(setdiff(unique(zoop_data$Station), unique(zoop_data_continuous$Station)))==6
## [1] TRUE
str(zoop_data_continuous)
## tibble [554,748 × 33] (S3: tbl_df/tbl/data.frame)
##  $ Source      : chr [1:554748] "EMP" "EMP" "EMP" "EMP" ...
##  $ SizeClass   : chr [1:554748] "Macro" "Macro" "Macro" "Macro" ...
##  $ Volume      : num [1:554748] 2.16 2.16 2.16 2.16 2.16 2.16 2.16 2.16 3.46 3.46 ...
##  $ Lifestage   : chr [1:554748] "Adult" "Adult" "Adult" "Adult" ...
##  $ Taxname     : chr [1:554748] "Alienacanthomysis macropsis" "Deltamysis holmquistae" "Hyperacanthomysis longirostris" "Mysida_UnID" ...
##  $ Taxlifestage: chr [1:554748] "Alienacanthomysis macropsis Adult" "Deltamysis holmquistae Adult" "Hyperacanthomysis longirostris Adult" "Mysida_UnID Adult" ...
##  $ SampleID    : chr [1:554748] "EMP NZ092 2015-03-12" "EMP NZ092 2015-03-12" "EMP NZ092 2015-03-12" "EMP NZ092 2015-03-12" ...
##  $ CPUE        : num [1:554748] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Phylum      : chr [1:554748] "Arthropoda" "Arthropoda" "Arthropoda" "Arthropoda" ...
##  $ Class       : chr [1:554748] "Malacostraca" "Malacostraca" "Malacostraca" "Malacostraca" ...
##  $ Order       : chr [1:554748] "Mysida" "Mysida" "Mysida" "Mysida" ...
##  $ Family      : chr [1:554748] "Mysidae" "Mysidae" "Mysidae" NA ...
##  $ Genus       : chr [1:554748] "Alienacanthomysis" "Deltamysis" "Hyperacanthomysis" NA ...
##  $ Species     : chr [1:554748] "macropsis" "holmquistae" "longirostris" NA ...
##  $ Undersampled: logi [1:554748] FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ Date        : POSIXct[1:554748], format: "2015-03-12" "2015-03-12" ...
##  $ Station     : chr [1:554748] "NZ092" "NZ092" "NZ092" "NZ092" ...
##  $ Chl         : num [1:554748] 3.18 3.18 3.18 3.18 3.18 3.18 3.18 3.18 33 33 ...
##  $ Secchi      : num [1:554748] 104 104 104 104 104 104 104 104 43 43 ...
##  $ Temperature : num [1:554748] 15.7 15.7 15.7 15.7 15.7 15.7 15.7 15.7 20.7 20.7 ...
##  $ BottomDepth : num [1:554748] 9.26 9.26 9.26 9.26 9.26 ...
##  $ Tide        : chr [1:554748] "High slack" "High slack" "High slack" "High slack" ...
##  $ Year        : num [1:554748] 2015 2015 2015 2015 2015 ...
##  $ Datetime    : POSIXct[1:554748], format: "2015-03-12 11:40:00" "2015-03-12 11:40:00" ...
##  $ Turbidity   : num [1:554748] NA NA NA NA NA NA NA NA NA NA ...
##  $ Microcystis : chr [1:554748] NA NA NA NA ...
##  $ pH          : num [1:554748] NA NA NA NA NA NA NA NA NA NA ...
##  $ DO          : num [1:554748] NA NA NA NA NA NA NA NA NA NA ...
##  $ SalSurf     : num [1:554748] 0.376 0.376 0.376 0.376 0.376 ...
##  $ SalBott     : num [1:554748] 0.404 0.404 0.404 0.404 0.404 ...
##  $ Latitude    : num [1:554748] 38 38 38 38 38 ...
##  $ Longitude   : num [1:554748] -121 -121 -121 -121 -121 ...
##  $ Month       : num [1:554748] 3 3 3 3 3 3 3 3 5 5 ...

Assess sampling density across space and time

zoop_surveys <- zoop_data_continuous %>% 
  select(Month, Year, Station, Longitude, Latitude) %>% 
  distinct(Month, Year, Station, Longitude, Latitude) %>% 
  mutate(Sampling=1)

zoop_emp <- zoop_surveys %>% 
  group_by(Station, Year, Longitude, Latitude) %>% 
  summarize(Sampling=sum(Sampling))

ggplot(zoop_emp, aes(Longitude, Latitude))+
  geom_point(aes(color=Station, size=Sampling), pch=21)+
  theme_bw()+
  facet_wrap(Year~.)+
  scale_color_viridis_d()

Fish + Zoop

Layer them on top of each other

# map for each year
ggplot()+
  geom_point(data=Fallmwt, aes(Longitude, Latitude), pch=21, color="gray50")+
  geom_point(data=Baystudy, aes(Longitude, Latitude), pch=21, color="orange")+
  geom_point(data=zoop_emp, aes(Longitude, Latitude), pch=4, color="slateblue")+  
  theme_bw()+
  facet_wrap(Year~.)+
  labs(subtitle="Spatial overlap - annual")

# map for all years
ggplot()+
  geom_point(data=Fallmwt, aes(Longitude, Latitude), pch=21, color="gray50", size=3)+
  geom_point(data=Baystudy, aes(Longitude, Latitude), pch=21, color="orange", size=3)+
  geom_point(data=zoop_emp, aes(Longitude, Latitude), pch=3, color="slateblue", size=3)+  
  theme_bw()+
  labs(subtitle="Spatial overlap - overall")

Water quality

Same stations as EMP zooplankton that are already plotted. Data lives here.

# library(discretewq)
# wq_stations <- wq(Sources = "EMP")
# wq_stations <- pivot_longer(wq_stations, 
#                             cols = c(Temperature, Chlorophyll:TKN, Salinity), 
#                             names_to = 'Analyte', values_to = 'Value')
# detach(package:discretewq, unload=TRUE)

Benthic

Code borrowed from the benthic.Rmd by JC

ben <- read_csv("data_raw/DWR Benthic raw count data 1975-2020 2021_10_01.csv")  %>% 
  mutate(StationClean =  gsub("-[A-Z]","",StationCode)) %>% 
  mutate(mon = lubridate::ym(paste0(Year, Month, sep = "-"))) %>% 
  group_by(StationClean) %>% 
  mutate(Latitude = mean(Latitude, na.rm = T), 
         Longitude = mean(Longitude, na.rm = T)) %>% 
  ungroup()
sites_mon <- ben  %>% 
  distinct(StationClean, Latitude, Longitude, mon) 

sites <- ben %>% 
  group_by(StationClean, Latitude, Longitude) %>% 
  summarise(n_months = length(unique(mon)), .groups  =  "drop") %>% 
  st_as_sf(coords = c('Longitude','Latitude'), crs = 4326,  remove = FALSE)

Sites that were sampled each month

Table of site information

Core Stations

core <- c("P8",
          "MD10",
          "D8",
          "D7",
          "D6",
          "D41",
          "D4",
          "D28A",
          "D26",
          "C3",
          "C10")

core_data <-  sites_mon %>% 
  filter(StationClean %in% core)

sites_core <- sites %>% 
  filter(StationClean %in% core)

Pesticide data

# library(dataRetrieval)
# ?dataRetrieval

Annual map

Regional designations

library(deltamapr)